package org.cosmocalc;

public class Trapezoid {
	public static double trapzd(MathFunction func, double a, double b, int n) {
		if (n == 1) {
			s = 0.5 * (b - a) * (func.f(a) + func.f(b));
			return s;
		} else {
			int it = 1;
			for (int j = 0; j < n - 2; j++)
				it *= 2; // Subdivisions
			double tnm = it;
			double delta = (b - a) / tnm; // Spacing of points to be added
			double x = a + 0.5 * delta;
			double sum = 0.0;
			for (int j = 0; j < it; j++) {
				sum += func.f(x);
				x += delta;
			}
			s = 0.5 * (s + (b - a) * sum / tnm);
			return s;
		}
	}

	private static double s; // Holds current estimate of integral. Do only 1
								// at a time.
}